      SUBROUTINE XFIND(ROD,VOD,ZMU,TST,TIME,XSAT,XSATDT)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/XKPLER/AS,E,F,EI,BW,W,BWDOT,WDOT
C
      COMMON/CFINDX/ AX2CHK,TCHECK,LCHECK
C
C                'XFIND' IS THE INTERNAL (KEPLERIAN) ORBIT GENERATOR.
C                XFIND DETERMINES THE SATELLITE POSITION AND
C                VELOCITY AT INPUT TIME.
C
C
C           INPUT
C                 ROD  SATELLITE POSITION AT START TIME
C                 VOD  SATELLITE VELOCITY AT START TIME
C
C           OUTPUT
C                 XSAT    SATELLITE POSITION AT SIMULATION TIME
C                 XSATDT  SATELLITE VELOCITY AT SIMULATION TIME
C
C           NO SUBROUTINES
C           NO TAPES
C
C
      DIMENSION R(3),VD(3),BR(3,3),BV(3,3),TM2(3),TM3(3),ROD(3),VOD(3),
     .          XSAT(3),XSATDT(3)
      DIMENSION RA(3,3),RST(3,3),RSULT(3,3),DUMY(3,3)
      DIMENSION RRR(3),VVV(3)
C
C
C
      T24=24.D0
      T72=720.D0
      T40=40320.D0
      T16=1.D0/6.D0
      T12=120.D0
      T50=5040.D0
      T36=362880.D0
C
      RODTVO=0.D0
      ROMD=0.D0
      VOSQD=0.D0
      DO 10 I=1,3
      VOSQD=VOSQD+VOD(I)**2
      ROMD=ROMD+ROD(I)**2
   10 RODTVO=RODTVO+ROD(I)*VOD(I)
C
      ROMD=DSQRT(ROMD)
C
      ALFAD=2.D0/ROMD-VOSQD/ZMU
      SQMU=DSQRT(ZMU)
      X=XPAST
   20 L=1
      TEST=TIME-TST
      IF(TEST.NE.0.D0) GO TO 30
      DN=0.D0
      X=0.D0
      DTDX=ROMD/SQMU
      X2CFN=0.D0
      X3SFN=0.D0
      GO TO 90
   30 X=X+(TEST-DPAST)/DTDX
      GO TO 50
   40 DTDX=RODTVO*(X-ALFAD*X3SFN)/ZMU
     .    + ((1.D0-ROMD*ALFAD)*X2CFN + ROMD)/SQMU
      X=X-(DN-TEST)/DTDX
   50 ALFX2=ALFAD*X**2
      DIFF=DABS(ALFX2)-AX2CHK
      IF(DIFF.GE.0.D0) GO TO 60
      XSQR=X**2
      XCUBE=XSQR*X
      ALFX22=ALFX2*ALFX2
      ALFX23=ALFX2*ALFX22
      X2CFN=(0.5D0 - ALFX2/T24 + ALFX22/T72 - ALFX23/T40)*XSQR
      X3SFN=(T16 - ALFX2/T12 + ALFX22/T50 - ALFX23/T36)*XCUBE
      GO TO 80
C  60 IF(SNGL(ALFAD).GE.0.) GO TO 70
 60   IF(ALFAD.GE.0.0D0) GO TO 70
      ARG=DSQRT(-ALFX2)
      EXPLS=DEXP(ARG)
      EXMIN=DEXP(-ARG)
      X2CFN=-(0.5D0*(EXPLS+EXMIN)-1.D0)/ALFAD
      ALFFF=DABS(ALFAD)
      X3SFN=(0.5D0*(EXPLS-EXMIN) - ARG)*((ALFFF)**(-1.5))
      GO TO 80
   70 ARG=DSQRT(ALFX2)
      X3SFN=(ARG-DSIN(ARG))*ALFAD**(-1.5D0)
      X2CFN=(1.D0-DCOS(ARG))/ALFAD
   80 DN=RODTVO*X2CFN/ZMU + ((1.D0-ROMD*ALFAD)*X3SFN + ROMD*X)/SQMU
      DIFF=DABS(TEST-DN)-TCHECK*TEST
      IF(DIFF.LE.0.0D0) GO TO 90
      L=L+1
      IF(L-LCHECK.LE.0) GO TO 40
      IEXIT=2
      RETURN
   90 DPAST=TEST
      XPAST=X
      RM=0.D0
      DO 100 I=1,3
      R(I)=(1.D0-X2CFN/ROMD)*ROD(I)+(TEST-X3SFN/SQMU)*VOD(I)
  100 RM=RM+R(I)**2
C
      RM=DSQRT(RM)
C
      DO 110 I=1,3
  110 VD(I)=SQMU*(ALFAD*X3SFN - X)*ROD(I)/(RM*ROMD)
     .     + (1.D0 - X2CFN/RM)*VOD(I)
C
      TDAY=TEST/86400.0D0
      BWST=BW+BWDOT*TDAY
      WST=W+WDOT*TDAY
      CALL TRANS(BW,EI,W,RA)
      CALL TRANS(BWST,EI,WST,RST)
      CALL MPYMAT(RST,RA,DUMY,1,2,RSULT,DUMY)
      CALL DMATVC(RSULT,R,RRR)
      CALL DMATVC(RSULT,VD,VVV)
      DO 120 I=1,3
      XSAT(I)=RRR(I)
  120 XSATDT(I)=VVV(I)
C
C   END OF MODIFICATION
C
C
C
      RETURN
      END
